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3D weak lensing with spin wavelets on the ball 
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We construct the spin flaglet transform, a wavelet transform to analyze spin signals in three dimen¬ 
sions. Spin flaglets can probe signal content localized simultaneously in space and frequency and, 
moreover, are separable so that their angular and radial properties can be controlled independently. 
They are particularly suited to analyzing of cosmological observations such as the weak gravitational 
lensing of galaxies. Such observations have a unique 3D geometrical setting since they are natively 
made on the sky, have spin angular symmetries, and are extended in the radial direction by addi¬ 
tional distance or redshift information. Flaglets are constructed in the harmonic space defined by 
the Fourier-Laguerre transform, previously defined for scalar functions and extended here to signals 
with spin symmetries. Thanks to various sampling theorems, both the Fourier-Laguerre and flaglet 
transforms are theoretically exact when applied to bandlimited signals. In other words, in numerical 
computations the only loss of information is due to the finite representation of floating point num¬ 
bers. We develop a 3D framework relating the weak lensing power spectrum to covariances of flaglet 
coefficients. We suggest that the resulting novel flaglet weak lensing estimator offers a powerful al¬ 
ternative to common 2D and 3D approaches to accurately capture cosmological information. While 
standard weak lensing analyses focus on either real or harmonic space representations {i.e., cor¬ 
relation functions or Fourier-Bessel power spectra, respectively), a wavelet approach inherits the 
advantages of both techniques, where both complicated sky coverage and uncertainties associated 
with the physical modeling of small scales can be handled effectively. Our codes to compute the 
Fourier-Laguerre and flaglet transforms are made publicly available. 

PACS numbers: 95.36.+x, 98.62.Py, 98.80.Es, 02.70.Hm, 02.70.Rr 


I. INTRODUCTION 

Weak gravitational lensing by large-scale structure, for example cosmic shear, is an observational probe that has 
the potential to constrain both the geometry and growth of structure of the Universe. It is a sensitive probe of 
dark energy physics, extensions to general relativity, and neutrino mass and hierarchy (see Refs. [1-5] and references 
therein for more details). The basic measurements of weak gravitational lensing from galaxy surveys are the third- 
flattening, or third-eccentricity, of galaxy images (colloquially refereed to as ‘ellipticity’), and galaxy sizes. These 
contain information about the intrinsic unlensed shape of the galaxies and the additional ellipticity (called ‘shear’), as 
well as size changes caused by the weak gravitational lensing effect along the line of sight. This angular information 
can be supplemented with the redshift of the galaxies (either photometric or spectroscopic), yielding a catalogue of 
galaxy positions, shapes, sizes and redshifts. These can then be compared to predictions from cosmological models in 
order to constrain their parameters. 

Due to the cosmological principle, i.e., assumptions of isotropy and homogeneity, the mean of the cosmic shear 
averaged over all galaxies in a sufficiently large survey is expected to be zero. However, the power spectrum of the 
shear (pair correlations of modes in 3D Fourier space) is expected to be nonzero and depends on both the geometry 
of the Universe and the power spectrum of matter perturbations in the 3D volume within which the galaxy sample 
lies. The 3D power spectrum of the galaxy shear estimates themselves is known as ‘3D cosmic shear’. There exist 
various strategies and approximations to efficiently compute the 3D cosmic shear power spectrum, including linking 
angular and radial scales, and binning in redshift—known as ‘tomography’—that are more widely used than the full 
3D case because of computational ease. 

From a theoretical perspective, 3D cosmic shear has been developed in a series of papers [6-9] that introduced the 
relation between the underlying matter power spectrum, the geometric lensing kernel, and the on-sky measurements 
of galaxies shear and size changes. Because the gravitational lensing signal is caused by a tidal effect around massive 
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objects, the amount of deflection that a photon experiences is related to the second derivative of the local gravitational 
potential, integrated along the line of sight. The resulting 3D shear distribution is thus well characterised by a 3D 
cosmic shear power spectrum defined in Fourier-Bessel space (the harmonic space defined by the eigenfunctions of the 
spherical Laplacian, see also [10] for a more generic presentation). While the theoretical side of 3D cosmic shear is 
well established, the application to data and inference of cosmological parameters are less developed; in fact, the 3D 
cosmic shear framework has only been applied to data twice. In Ref. [11] a 3D cosmic shear analysis was performed 
using data observed on a small field of approximately one square degree area of sky, with the aim of measuring erg (the 
amplitude of matter fluctuations on 8 Mpe scales) and Om (the dimensionless matter density) as a proof of concept. 
In Ref. [12], a similar analysis was performed on the CFHTLenS survey [13] and used to constrain cosmological 
parameters including the dark energy equation of state, which is parametrised as w{z) = wq ^ Waz/{1 + z). On 
large scales, > 1.5h“^Mpc, it was found that results were consistent with the measurements of the cosmic microwave 
background (CMB) anisotropies from the Planck satellite [14]. But on small scales, < 1.5h“^Mpc, results were found 
to be inconsistent and favoured less clustering of matter than that predicted by Planck, with the interpretation that 
this could be due to baryonic feedback effects or systematics in the data. 

Due to various observational constraints, galaxy surveys cover small, and potentially disconnected, regions of the 
sky. This leads to the problem that the observed power spectrum is related to that of the full sky through a 
convolution with the window function of the survey. Furthermore, inhomogeneous observation strategies, and galaxy 
sample selection can result in complicated three dimensional masks or weight maps, which significantly increase the 
computational requirements of 3D methods. In fact, even 2D masks will mix angular and radial modes in a nontrivial 
way [12]. Two methods have been proposed to take such masks into account. The first is a pseudo-C^ approach where 
the effect of the mask on the Fourier-Bessel modes is computed in a mixing matrix. Then, either the theory is forward 
convolved, or the data is deconvolved with the inverse mixing matrix. This has been applied to data only once in the 
weak lensing context but is well known from power spectrum measurements of the CMB and galaxy clustering (see, 
e.g., Refs. [15-18]). Such a process is computationally difficult as the mixing matrices tend to be close to singular - for 
weak gravitational lensing surveys the mask is also very complex on small scales due to masking of stars, cosmic rays 
and other image artifacts. The second approach is through a forward modeling of the data using Bayesian hierarchical 
modeling [19], where the mask is included in the model and the data given infinite variance in the regions of the mask 
(similar to the approach used to handle a mask by [20]). This approach provides the full posterior probability of 
model parameters but is computationally expensive. 

Another area of difficulty in the analysis of the cosmic shear power spectrum is the interpretation of scale-dependent 
features in the power spectrum. For a given power spectrum of matter perturbations the cosmic shear power spectrum 
is readily computable from theory. However, because of nonlinear structure growth {e.g., Ref. [21]), baryonic feedback 
effects {e.g., Ref. [22]), photometric redshift systematics, and uncertainty in modeling neutrino physics {e.g.. Ref. [23]), 
the modeling of the matter power spectrum on small-scales of less than ^ lh“^Mpc is highly uncertain. Therefore, 
designing statistics and weighting of the data that mitigate the use of the uncertain small-scale regime is particularly 
important in extracting parameter estimates robustly. 

In this paper we deal with the full 3D cosmic shear case. However, it should be noted that the majority of theoretical 
papers on weak gravitational lensing power spectra use cosmic shear ‘tomography’, and all but two papers [11, 12] 
that use cosmic shear for cosmological parameter estimation do not use power spectra but instead use real-space 
correlation function measurements. Using correlation functions from the data is ostensibly a good way of accounting 
for the mask in the data because the measured correlation function is weakly affected by it. However, the covariance 
properties of the data are much more complicated {e.g.. Refs. [24-26]) and result in the need to handle the survey 
mask accurately due to ‘beat coupling’. Furthermore, the kernel through which correlation functions depend on the 
power spectrum of matter perturbations is different for each real-space angular measurement, for each redshift bin, 
and can extend to very small-scales {e.g.. Ref. [27]). In this paper we explain how these issues may be alleviated by 
working in wavelet space. 

Wavelets are basis functions or kernels that are well localized in both real and frequency space. They are used in a 
number of disciplines to solve challenging data analysis, data compression, feature extraction, and pattern recognition 
problems. In cosmology, wavelets on the sphere have become an integral part of state-of-the-art algorithms to analyze 
data from the CMB, in particular to isolate the CMB signal from astrophysical foregrounds, extract its statistical 
properties, and search for potential anomalies (see, e.g.. Refs. [28-37] and references therein). As with the CMB, 
observations of galaxy surveys are made on the celestial sphere. Yet, wavelets on the sphere are not typically used in 
standard analyzes of galaxy surveys because the cosmological information of interest is in the full 3D distribution of 
galaxy properties; ignoring the radial information causes a significant loss of sensitivity to the geometry and growth 
of structure of the Universe. Various wavelet methods have been developed to analyze galaxy survey data recently 
but none of them are suited to 3D weak lensing observables, or 3D signals with spin angular symmetries in general. 
As mentioned previously, weak lensing observables such as cosmic shear are spin signals on the sky, extended to three 
dimensions through additional radial information. Thus, uncertainties and observational complications are separable 
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on the sky and along the line of sight {e.g., the ability to distinguish galaxies from stars and to estimate their redshifts). 

A number of 2D and 3D methods, including 3D wavelets [38-40], have been developed to analyze projected weak 
lensing observables as well as the polarization of the CMB, which is also a spin signal {e.g., Refs. [41-45]). However, 
these approaches do not natively deal with 3D spherical geometry and spin symmetries simultaneously. Wavelets in 3D 
spherical coordinates were constructed in Ref. [38] using B-spline kernels and the Fourier-Bessel basis. While these are 
directly connected to the Fourier-Bessel formalism, they can probe only 3D isotropic features in scalar signals. Hence, 
they do not account for the angular spin symmetries or the unique 2+ID geometry of weak lensing observables. By 
contrast. Ref. [39] introduced ‘flaglets’ to probe the angular and radial scales separately. This approach was based on 
the Fourier-Laguerre transform, a novel separable 3D harmonic transform with a sampling theorem and an analytical 
connection to the Fourier-Bessel basis. More recently. Ref. [40] proposed a construction of separable 3D needlets 
based on the Fourier transform in the radial direction. 

In this paper, we extend the framework of Ref. [39] and present novel Fourier-Laguerre and flaglet transforms 
supporting spin symmetries and also directional features in the angular direction. We detail sampling theorems 
and efficient algorithms to compute these transforms exactly, exploiting the recently constructed spin directional 
wavelet transform on the sphere and sampling theorem on the rotation group, and their corresponding fast algo¬ 
rithms, developed in Refs. [46-48]. Thanks to the separability of the radial and angular components, the novel spin 
Fourier-Laguerre and flaglet transforms natively support the 2+ID spherical geometry of the weak lensing data while 
accounting for angular spin symmetries. They can also be related to the standard 3D weak lensing framework in the 
Fourier-Bessel basis thanks to the analytical connection between the Fourier-Laguerre and Fourier-Bessel transforms. 
Finally, we discuss the advantages of working in flaglet space (instead of the standard correlation functions, Fourier, 
or Fourier-Bessel approaches) that arise from the dual localization properties of flaglets in pixel and frequency space, 
for example in dealing with complicated angular masks or small scales where physical modeling is less certain. 

This paper is structured as follows. In Sec. H we give a short review of the 3D weak lensing formalism. We adopt 
the language of observational cosmology and focus on the observables that can be measured from galaxy survey data 
and connected to cosmological models. In Secs. HI and IV we step back from this formalism and give a formal 
mathematical presentation of the new analysis techniques (transforms and operators) for radial, spherical and 3D 
fields. This leads to the definition of the novel spin Fourier-Laguerre and spin flaglet transforms (in Secs. HI and IV, 
respectively). In Sec. V, we combine the two viewpoints, apply the spin flaglet transform to 3D weak lensing, and 
study the properties of the shear flaglet coefficients. We present conclusions in Sec. VI. Further technical details on 
the special functions and useful approximations used in this paper are presented in appendixes. 


II. 3D WEAK LENSING 

In this section we give a brief introduction to the 3D weak lensing formalism, which has been developed in a 
series of recent papers [6-9]. Here we focus on the basic principles of 3D weak lensing and discuss the lensing 
potential, observable quantities, and how observables and theory can be compared to constrain cosmological models 
and parameters. For further background we refer the reader to the excellent exposition of 3D weak lensing given in 
Ref. [7]. Details of practical aspects and data-related complications are addressed in Ref. [12] {e.g., pseudo-C^ methods, 
real and imaginary covariance structures, E and B-mode decomposition, and weighting due to shape measurement 
biases). We touch these issues briefly in this section, and in Sec. V, but leave detailed investigations on data and 
simulations to future work. 


A. Lensing potential 


The weak gravitational lensing effect is commonly expressed in terms of the lensing potential 0, which depends on 
the integrated deflection angle along the line of sight, sourced by the local Newtonian potential 4>, 


4>{r,n) = 


f 


dr 


, fKjr-r') 
fK{r)fK{r') 




( 1 ) 
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where r is the comoving distance, n the angular position on the sky, and c is the speed of light in a vacuum. The 
geometrical factor reads 


fK{r) 


' sin(r), if K = 1 

< r, if K = 0 

, sinh(r), if K = —1 


( 2 ) 


for cosmologies with positive, flat and negative global curvatures, with curvature = 1, 0, and —1, respectively. This 
expression assumes the linear regime and also the Born approximation, i.e., that the path of photons is unperturbed 
by lenses. 

The Newtonian potential 4> can be related to perturbations in the underlying matter density 6 via Poisson’s equation. 


V^4>(r, n) 


2a(r) 


S{r, n), 


(3) 


where flu is the dimensionless matter density, Hq is the current value of the Hubble parameter, and a(r) is the 
dimensionless scale factor. The 3D gradient is defined relative to comoving coordinates. This relation can be expressed 
conveniently in terms of the Fourier-Bessel basis since the Fourier-Bessel basis functions are eigenfunctions of the 
Laplacian operator V^, with eigenvalues —k^. Consequently, the harmonic representation of Eq. (3) reads 

^em{k-,r) = (4) 

where £ and m label angular harmonic modes and k labels the radial wavenumber (in units of h~^ Mpc). The 
Fourier-Bessel transform is defined expiicitiy in Sec. HID by Eqs. (37) and (38), and discussed extensiveiy in the same 
section. The dependency on r must be retained in Eq. (4) to account for the time evoiution of the fieid: we appiy 
the Eourier-Bessei transform to the homogenous fieid existing everywhere at the cosmoiogicai time corresponding to 
comoving distance r. 

The harmonic representation of the lensing potentiai then reads 

<Pem{k) = ^ f drr'^ji{kr) f f (5) 

7rc2 7ji+ Jo fK{r)fK[rJ Jr+ 


where j£ are the sphericai Bessei functions and M+ denotes the positive reai haif-iine, i.e., [0,oo). The difference 
between Eq. (5) and the expression shown in Ref. [7] is due only to the different conventions used for the spherical 
Bessel transform, where we adopt the convention of Refs. [6, 39] (see Appendix A and Sec. HID). 


B. Observables 


Weak lensing generates distortions in the observations of a background field, which may be characterised by spin 
quantities of several orders (see, e.g., Ref. [49]), with spin number s. In the weak lensing regime (z.e., away from the 
critical curve of lensing masses, where there are not multiple images of sources) four distortions can be produced [49]: 
the size magnification qk [50]; the shear 27 ; and the flexion, the combined effect of the first-flexion iJ^ (a centroid 
shift) and third-flexion (a trefoil distortion). We introduce the generic notation to denote these spin quantities, 
which are related to the lensing potential 0 by (see, e.g., Ref. [7]) 


sX{n,r) 


iJ’(n,r) 

< 

2l{ri,r) 

< sGin.r) 


I(55 + 55)0(n,r), if s = 0 

^(555 + 555 + 555)0(n,r), if s = 1 
^(55)(/)(n,r), if s = 2 

^(555)(/)(n,r), if 5 = 3 


(6) 


where 5 and 5 are spin raising and spin lowering operators [51-54]. Eurther details about these operators and spin 
spherical harmonics in general are given in Sec. Ill and Appendix A. We now focus on the spin -2 quantity of shear 
27 , which is of considerable interest since among the spin quantities of it has the highest signal-to-noise ratio in 
observational data (see, e.g., Refs. [55, 56]). We revisit the other spin quantities of ^y in Sec. VB (where we present 
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their flaglet covariances). 

The shear field ±27 is typically decomposed into its real and imaginary parts, ±27 = 7i ± 172 , and is related to the 
leasing potential by [7] 


27 (n, r) = 7 i (n,r)+i 72 (n,r) = + i^!'^(n,r))/2, (7) 

27*(n,r) = _27(n,?’) = 7 iK 0-i72(n,r) = - i?!'^(n,r))/2, (8) 

where * denotes complex conjugation and we have represented the leasing potential by its parity even and odd 
components given by the E- and B-mode signals, respectively, denoted by superscripts E and B. Alternatively, the 
E- and B-mode signals may be related to the shear by 


= {Q‘^2l{n,r)+Q‘^-2j{n,r)), (9) 

4>^in,r),= -i{S‘^ 2 'l{n,r) - d‘^- 2 'y{n,r)), (10) 

where and are normalized versions of and and, when expressed in their Eourier-Bessel representations 
(see Sec. HID), are related by 


= (7V,,2)7fn(fc), (11) 

= (iv,,2)7L(fc), ( 12 ) 

where N (^^2 = \J [72)! ~ (-^^,- 2 )”^- The shear induced by gravitational lensing produces an E-mode signal only since 
density (scalar) perturbations cannot induce a parity odd B-mode component [7]. In the absence of systematic effects 
(see, e.g., Refs. [55, 56]), = 0 and = 0. The 3D cosmic shear formalism nevertheless includes B-modes since 

the B-mode signal—and EB-mode cross correlations—computed from data can be used as a null-test to search for 
residual systematics. 

The spin ±2 cosmic shear signal ±27 can also be represented in Eourier-Bessel space, which leads to a simple 
harmonic connection between the shear and the lensing potential. Since shear is a spin quantity, it is decomposed 
into its spin Eourier-Bessel coefficients (see Sec. HID), denoted ± 27 ^m(^)- Scalar harmonic E- and B-mode signals 
may be computed from ± 27 ^m(^) by 


Eem{k) = - {2jem{k) + -27em{k))/2, (13) 

Blm{k) = i(27«m(fc) - -2llm{k))l2, (14) 

where yet another normalization is assumed (as is standard practice). It also follows that 

±27^m(^) = - ±i^^m(^))- (16) 

By Eqs. (13) and (14) and either Eqs. (7) and (8) or Eqs. (9) and (10), it follows that the lensing potential and the 
shear are related in harmonic space indirectly by 

<^f„(fc) = -2Ne,-2EeM, (16) 

4>f^ik) = -2Ne,_2Bem{k). (17) 

Consequently, the shear can be related directly to the lensing potential in Eourier-Bessel space by 

27«m(fc) = ^Ne,24>fm{k)- (18) 

In fact, Eq. (18) can be recovered directly from the harmonic representation of Eq. (6). Nevertheless, we expose the 
E- and B-mode decomposition since it is the standard approach and is of additional practical use in studying residual 
systematics, as outlined. 

Einally, note that galaxy distances are not directly measurable; only their redshifts can be estimated from galaxy 
colours or spectra. Thus, the 3D shear signal ± 27 (^) in real space is typically constructed from redshift space 
observations via the relation r{z), which assumes a fiducial reference cosmological model [12]. 
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C. Cosmology 

Eqs. (4) and (5) show that the leasing potential contains information about the matter fluctuations and geometry of 
the Universe, via S and Fk^ respectively. In the former, most of the information arises from the first nonzero moment 
of the field, the 2-point covariance of the the fluctuations, described by the matter power spectrum P(/c;r): 

( {k - k'), (19) 

where the angle-brackets (without a bar) indicate an ensemble average over several realizations of the field given the 
same cosmological model and power spectrum P(/c;r), denotes the ID Dirac delta function and 6^ the Kronecker 
delta symbol. Note that we have assumed homogeneity and isotropy and that the field is evaluated at a fixed comoving 
distance [7]. Higher order statistics of the density field are of less interest because they are more difficult to model 
and measure in data. For this reason we focus on the 2-point power spectrum of the cosmic shear and other weak 
leasing observables. 

The 3D leasing power spectrum Cf^{k^ k') is defined by 

( 4>fm{k) ik') ) = cf (fc, k')5Ul^,, (20) 

and is dependent on P{k;r) and P^- Note that the 3D leasing potential is not homogenous and isotropic since it is 
given by a 2D projection of the Newtonian potential between an observer and source at radius r. Consequently, the 
leasing potential is homogenous and isotropic in the angular direction but not in the radial direction [7] . 

By Eq. (18), and equivalently by Eq. (6) for s = 2, the covariance of the 3D cosmic shear 27 in Fourier-Bessel space 
is related to the leasing power spectrum by 

{2lUk)2ll^,{k')) = \{N,,2?Ct^{k,k')5l,5l^,. ( 21 ) 

Note that since the shear transform coefficients are complex quantities, their joint probability distribution must 
include the correlation between the real and imaginary parts [57, 58]. This may be accounted for by creating an 
‘Affix’ covariance, as shown in Ref. [12], or by adopting real spherical harmonics. 

Further observational aspects such as the inclusion of photometric redshift information or a varying number-density 
of sources can also be included as shown in Ref. [9]. A further complicating aspect is that the observed galaxy 
ellipticity is not a direct measure of the shear but is a combination of the shear and the intrinsic (unlensed) ellipticity 
of the galaxy: e°^®(n,r) = e^(n,r) + 7(71, r), to linear order. These effects can be taken into account by modeling 
auto-correlations between the intrinsic ellipticities e^(n, r) and cross-correlations between the intrinsic ellipticities and 
the shear as described in Ref. [59] (correlations with CMB lensing can also be included, as also shown in Ref. [59]). 
In this paper we are concerned with linking observable power spectra to the model power spectrum that should be 
assumed to include all of these effects. However, for simplicity and clarity we will refer to only the ‘shear’ power 
spectrum and keep the notation Cf^{k^ k'). For more details on the additional modeling aspects we refer the reader 
to the references provided. 

We do not consider higher-order statistics, such as the bi-spectrum or the tri-spectrum [59]. These can be expressed 
as generalizations of the power spectrum method presented. Note that in the Fourier-Bessel transforms we have 
assumed a flat geometry, and that the more general case requires the use of ultra spherical Bessel functions, where the 
geometry is now a generalized manifold. However, as pointed out in Ref. [11], the use of the spherical Bessel functions 
is well motivated in the limit of 1 and k ^ (curvature scale) 

The standard approach of 3D cosmic shear, as performed in Ref. [12], is to measure the Fourier-Bessel coefficients 
of the shear 2 lem{k) from a galaxy survey and then use them to constrain Cf^{k^ k') (hence, P{k] r) and Fk) through 
Eq. (21). This harmonic-space approach can be challenging due to both observational and technical complications. As 
mentioned previously, it is difficult to simultaneously deal with complicated masks, mode coupling and scale mixing 
of the Fourier-Bessel coefficients. We will come back to these issues in Sec. V when we detail our new method to 
constrain from flaglet coefficients. 

We now step back from weak lensing and adopt a more general signal processing viewpoint to define the novel 
3D spin Fourier-Laguerre and flaglet transforms (which are suitable for the analysis of arbitrary 3D spin fields). In 
Sec. V we apply these transforms specifically to weak lensing observables and construct an estimator related to the 
3D lensing power spectrum. 
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III. SPIN FOURIER-LAGUERRE TRANSFORM 

After concisely reviewing signals and transforms defined on the radial line IR+, sphere and rotation group 
SO(3), we construct an exact, separable Fourier-Laguerre transform to analyze signals of arbitrary spin defined on 
§^x M+, i.e., the 3D space formed by the product of the sphere and the radial line M+. To subsequently construct a 
directional wavelet transform in this space (presented in Sec. IV), it is necessary to also consider signals and transforms 
defined on SO( 3 )x]R+, i.e., the space formed by the production of the rotation group SO(3) and the radial line IR+. 
We thus also construct a Wigner-Laguerre harmonic transform on SO(3) xM+. Sampling theorems and fast algorithms 
to compute these transform exactly for bandlimited signals are presented. In addition, the translation and convolution 
operators needed to construct a 3D spin wavelet transform are outlined. Finally, we derive the relation between the 
Fourier-Laguerre and Fourier-Bessel transforms, allowing one to exactly compute the Fourier-Bessel coefficients of 
signals bandlimited in Fourier-Laguerre space. This property is of use in Sec. V to connect the Fourier-Laguerre and 
flaglet transforms with the 3D weak leasing formalism. Note that in what follows brackets with a separating bar 
denote inner products (and not ensemble averages). In other words, {f\g) = /^d^;/^*, where the integral runs over 
the space S of interest, i.e., where / and g are defined, with measure d^;. 


A. Transforms on R+, and SO(3) 


We first consider a square integrable complex signal on the radial line / G I/^(M+), which can be expanded using 
the spherical Laguerre transform introduced in Ref. [39] , with forward and inverse. 


fp = (/l-f^^p) = [ 

Jr+ 

OO 

/(^) = 

p=0 

respectively. The Laguerre basis functions are specified by 


( 22 ) 

(23) 


Kp{r) = 


,-r/2T 


{p + a)\ A/r«+i 





(24) 


where is the p-th generalized Laguerre polynomial of order a (see Appendix A for further details about the 
basis functions Kp). In this paper we consider the case a = 0. Note that r G is a scaling factor to map the M+ 
sampling theorem to any finite interval [0, R] of interest. We refer the reader to [39] for more details on the spherical 
Laguerre transform. We take r^dr as the natural measure on M+ since we aim to construct 3D transforms in spherical 
coordinates, where the volume-invariant measure involves r^dr (in any case, an alternative measure could be adopted 
if it were desired). 

We now consider a square integrable complex function on the sphere s/(n) G I/^(§^) with n = (^,0), where 
0 G [0, tt] denotes colatitude and (j) G [0, 27r) longitude. A spin function sf with spin number s G Z transforms under a 
local rotation by 'd G [0, 27r) in the tangent plane centered on n as s/'(n) = e“'^^^s/(n), where the prime denotes the 
rotated function [51-54]. For s = ± 2 , such as the cosmic shear field, sf is invariant under local rotations of ± 7 r. Note 
that the sign convention adopted for the argument of the complex exponential differs from the original definition [52] 
but is identical to the convention used in the context of the polarization of the CMB [53, 54]. The natural harmonic 
transform on the sphere that accounts for spin symmetry is the spin spherical harmonics transform, with forward and 
inverse. 


sfim = isflsyim) = [ dn{n) sf{n) (25) 

OO £ 

sf{n) = Z] V sYimin), (26) 

£=0 m=—£ 

respectively, where denotes the spin spherical harmonics (using the Condon-Shortley phase convention; see 

Appendix A for more details). The usual invariant measure of the sphere is dfl{n) = sinOdOdcl). Note that spin 
signals have shm = O 5 < |s|. 

The spherical Laguerre and spherical harmonics transforms can be combined naturally into a 3D separable trans- 






FIG. 1: Real and imaginary parts of the Fourier-Laguerre basis function syim{'n)Kp{r) with s = 0, = 10, m = 5, p = 40, 

shown on the 3D Fourier-Laguerre sampling scheme with L = P = 80, plotted up to r = P/2 with P the radius of the final 
sample of the radial sampling scheme [39]. This plot shows that the 3D basis functions successfully probe harmonic modes 
in the angular and radial directions separately. Furthermore, notice that the basis functions are not well localized in space, 
similar to standard Fourier bases. 


form, as shown in Ref. [39] for the spin-0 case. Before we extend this construction to higher spin numbers in the next 
section, we turn to the rotation group, SO(3), and its standard harmonic transform: the Wigner transform [48, 60]. 
The rotation group is the natural manifold on which to construct a spherical wavelet transform probing directional 
features. This is due to the rotation properties of spherical harmonics and the natural convolution operator on the 
sphere (recalled in Sec. HIE). Thus, we consider a square integrable complex signal defined on the rotation group 
f{p) e L2(S0(3)), with p = ((a,;5,7), where a G [0,27r),/3 G [0,7r],7 G [0, 27r) are the Euler angles. Its forward and 
inverse Wigner transforms are given by 

fLn = {f\Dtn) = [ Mp) fiP) DinniP). (27) 

JS0(3) 

7 ^ 9 / 4 _ 1 _ ^ ^ 

/(/") = E E fLD^^nip), ( 28 ) 

£=0 m=—£n=—£ 

respectively, where are Wigner functions and the invariant measure on the rotation group reads d/i(p) = 

sin/3 dpdad^ (see Appendix A for further details). 


B. Fourier-Laguerre transform on and Wigner-Laguerre transform on SO(3)xR^ 

We now combine the previous transforms and extend the Eourier-Laguerre tranform presented in Ref. [39] so that 
its angular part can support spin signals on S^, in order to match the nature of weak leasing observations. More 
precisely, we now consider 3D complex functions in spherical coordinates s/(n, r) G I/^(S^xM+). We create a separable 
transform by combining the spherical Laguerre transform with the spin spherical harmonics, leading to the forward 
and inverse transforms 


sfem,p = {sf\ sYemKp ) = [ dn{n) [ drr\f{n,r) ,Yl^{n)K;{r), 

~ ^ ^ sf£m,p s^£m{'^) 

£mp 


respectively. The spin property implies that sf£m,p = 0, 


(29) 

(30) 
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For functions f{p,r) G I/^(SO(3)x]R+) we define the Wigner-Laguerre forward and inverse transforms as 

fL,p = if \ ) = [ Mp) [ drr^fip,r)Dl„{p)K;{r), (31) 

JS0(3) Jr+ 

fiP^^) = fUpDtnip) K,{r), (32) 

imnp 

respectively. In these expressions, the imn indices refer to the angular modes (in the spherical harmonics or Wigner 
transforms), while p corresponds to the radial mode (from the spherical Laguerre transform). An example of a basis 
function is shown in Fig. 1. Finally, as in the individual transforms, the summations for the harmonic modes run over 
p G N, G N, —< m < ^ and —i<n<i. In what follows, we will omit these bounds for conciseness and use the 
terminology ‘angular’ for both the and SO(3) parts of the signals and transforms. 


C. Sampling theorems and numerical implementation 

In practice, evaluating the integrals in Eqs. (29) and (31) numerically requires quadrature rules. While generic 
numerical integration methods can be adopted, this problem can be tackled more accurately by appealing to sampling 
theorems. When considering bandlimited signals, sampling theorems lead to exact quadrature rules, which allow one to 
discretize the space such that the previous integrals can be computed exactly. In practice, signals are transformed with 
accuracy close to the level of machine precision, with errors due to the (finite, approximate) representation of floating 
point numbers only. This approach is particularly desirable when performing numerous transforms successively or 
when high numerical accuracy is required (which is typically the case in modern observational cosmology). 

The radial and angular parts of the Fourier-Laguerre and Wigner-Laguerre transforms are fully separable. This 
implies that one can simply combine sampling theorems on the sphere, the rotation group and the radial line to 
obtain sampling theorems on §^x M+ and SO(3)xM+. We now recall these sampling theorems and the notation for 
the sampling points and quadrature weights. 

On M+, we use the spherical Laguerre quadrature rule presented in Ref. [39]. bandlimited signals have fp = 
0, Vp > P, with P the radial band-limit, while the sampling nodes and quadrature weights are denoted by and 
rc/c, respectively. This exact quadrature, with corresponding nodes and weights, is used to compute Eq. (22). 

On the sphere, we consider the sampling theorem developed in Ref. [61]: a function gf with band limit L satisfies 
sfim = 0, \/i > L and can be completely described by the values taken on the equiangular, separable pixelization: 
'f^ij = (^i: 0i) with i = 0,..., P — 1 and j = 0,..., 2P — 1. The spherical harmonics coefficients of Eq. (25) can be 
calculated with a discrete sum of sfi^^ij) with quadrature weights Wij. In principle, other sampling theorems could 
be used {e.g., Refs. [62, 63]), but note that at fixed band limit this sampling scheme requires the fewest samples [61]. 
pixelization schemes not based on sampling theorems could also be adopted {e.g., HEALPix [64]), in which case larger 
errors due to approximate numerical quadrature are expected. On SO(3), a function / with angular band limit L 
and azimuthal band limit N satisfies = 0, yi > P, |n| > N. The sampling theorem of Ref. [61] was extended 
to SO(3) in Ref. [48]: the sampling points, also equiangular and separable, are denoted by with 

i = 0,..., P — 1, j = 0,..., 2P — 1, and h = 0,..., 2N — 1. The quadrature weights are denoted by Whij and can be 
used to compute the Wigner coefficients of Eq. (27) exactly. 

The theoretically exact spin Eourier-Laguerre transform for bandlimited signals on S^x IR+ reads 

sfem,p = E sf{nij,rk) Kp{rk), (33) 

ijk 

sf{nij,rk) = E sfim,p sYemiriij) Kp{rk), (34) 

imp 


which are the discrete versions of Eqs. (29) and (30), respectively. 

The theoretically exact Wigner-Laguerre transform for bandlimited signals on SO(3 )xM+ reads 

fLn,p = E '^hij Wk f{phij,rk) Di,^{Pkij) Kp{rk), (35) 

hijk 

fiPhij,rk) = Y. fL,p D^r^niPhij) Kp{rk), (36) 

imnp 
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which are the discrete versions of Eqs. (31) and (32), respectively. 

Although we express the forward transforms of Eqs. (33) and (35) using discrete quadrature above, in practice we 
do not compute these expressions explicitly but rather compute them using fast algorithms, drawing on Refs. [48, 61], 
as discussed below. The sums in these equations are finite and adjusted to the band limits (for the imnp indices) 
and to the sampling nodes (for the hijk indices). In what follows we will omit the pixel indices and bounds of the 
summations for conciseness. We will also use the integral (continuous) forms of the transforms for clarity, since the 
sampling theorem guarantees that these can be evaluated exactly for bandlimited signals (which is typically the case 
in practical applications). 

As described above, any bandlimited signal with radial and angular band limits P and L, respectively, can be 
represented exactly by ^ 2PL^ samples on x IR+ thanks to the combination of the sampling theorem for the 
pixelization of Ref. [61] and the spherical Laguerre sampling theorem [39]. Similarly, signals on the rotation group, 
with the same band limits P and L, and an additional azimuthal band limit are captured by ^ APNL^ samples 
on SO(3)x]R+, thanks also to the corresponding sampling theorem on the rotation group [48]. 

In terms of computational complexity, the separability of the angular and radial components in the Eourier-Laguerre 
and Wigner-Laguerre transforms is an essential property. East algorithms exist to compute transforms on the sphere 
and the rotation group [47, 48, 61, 65-67]. In particular, we use the implementations detailed in Refs. [48, 61], which 
sets the complexity of the spin spherical harmonics and Wigner transforms to 0{L^) and 0{NL^), respectively, by 
exploiting fast Eourier transforms on the ring torus. The spherical Laguerre transform simply scales as 0(P), and the 
quadrature weights and basis functions are computed through recurrence relations [39]. Therefore, the spin Eourier- 
Laguerre and Wigner-Laguerre transforms scale as 0{PL^) and 0{NPL^)^ respectively. Note that the spin number 
5 is a parameter that does not change the complexity or the accuracy of any of the transforms, which is often not the 
case for alternative approaches. 


D. Connection to the Fourier-Bessel transform 

In contrast to the Eourier-Laguerre construction, the standard nonseparable basis for 3D spin functions s/(n,r) G 
I/^(§%<M+) (such as weak leasing observables) is the combination of spherical harmonics with spherical Bessel functions 
ji. These define the Eourier-Bessel transform, already used in Sec. II, with forward and inverse 

sfim{k) = { sf \ sYemje{k-) ) = [ dn{n)i - ( drr'^sf{n,r) sYf*^{n)j'^{kr), 

Jg2 V TT J^+ 

sf{n,r) = yZ\ - [ dkk"^ sfem{k) sYem{n)je{kr). 
im ^ ^ 

Unlike in the Eourier-Laguerre case, the Eourier-Bessel basis functions and transform are not separable: the angular 
modes i are coupled with the radial modes k through the spherical Bessel function ji(kr). This is because the basis 
functions are solutions of the (isotropic) Laplacian operator in spherical coordinates, with eigenvalues —It is 
important to note that the Eourier-Bessel transform is notoriously difficult to evaluate for generic functions because 
it does not admit a sampling theorem [39, 68]. Therefore, one must appeal to an approximate quadrature rule for the 
radial integrals (see e.g., Refs. [6, 38, 69]), which is challenging because the spherical Bessel functions have infinite 
support and oscillate rapidly. 

Exact analytical formula only exist for simple forms of signals, but the first exact formula to compute the Eourier- 
Bessel transform of a wide class of signals was presented in Ref. [39]. The starting point of this approach is to 
decompose the spherical Bessel functions in the spherical Laguerre basis, i.e., 


oo 


jeikr) = '^Je,p{k)Kp{r), 

p=0 

(39) 

Je,p{k) = f drr^Mkr)K;{r). 

Jr+ 

(40) 


(37) 

(38) 


The Eourier-Bessel coefficients of s/ can then be expressed in terms of its Eourier-Laguerre coefficients: 



V 


Je,p{k)- 


fe.m{k) = 


(41) 
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This sum is finite if gf is bandlimited in spherical Laguerre space, and sfim.p is also evaluated exactly thanks to the 
exact quadrature rules of the Fourier-Laguerre transform. Thus, Eq. (41) provides a way to compute the Fourier- 
Bessel transform of signals described in Fourier-Laguerre space exactly. It was also shown in Ref. [39] that 
admits an analytical form, involving the moments of ji{kr). Although Ji^p{k) can still be challenging to evaluate 
numerically for high it does not depend on the signal gf under consideration and can be tabulated. We will see 
that Ji^p{k) appears in the treatment of weak leasing observables, natively expressed in the Fourier-Bessel basis. 


E. Rotation, translation, and convolutions 

We now construct rotation and convolution operators on S^x IR+ and SO( 3 )x]R+, generalizing the operators of 
the scalar Fourier-Laguerre transform [39, 68 ]. In general, such operators are essential for constructing a meaningful 
wavelet transform, where the wavelet coefficients of a signal are its convolution with the wavelets. Since wavelets have 
compact support in real and frequency space, this approach allows the transform to extract well defined scales and 
directions. 

We first recall that the rotation of a spin function on the sphere s/ G is naturally defined by 

{TZp J){n) = (42) 

where x is the Cartesian vector corresponding to n and is the 3D rotation matrix corresponding to the rotation 
operator IZp. In harmonic space, the rotation conveniently reduces to [46] 

t 

(7^p = E <n(p) sftn- (43) 

n=—£ 

Furthermore, we define the translation of a radial function / G in spherical Laguerre space as [39, 68 ] 

{rrf)p = fpKp{r). (44) 

This is analogous to the harmonic representation of the translation operator for functions on the infinite (Cartesian) 
line M (but not analogous to the real space expression of the translation operator). In addition, such a translation 
operator can be expressed in terms of convolution with a shifted Dirac delta function, again in analog to translation 
on M, as detailed in Ref. [ 68 ]. 

We now consider functions on S^x M+, and introduce a directional convolution operator ® such that its action on 
/, ^ yields a function on SO( 3 )x]R+ defined by 

= if iTlpTt g) = [ dfJ(n) [ drr^ f{n,r) {TZpTt g)*{n,r). (45) 

Jr+ 

We also introduce an axisymmetric convolution operator 0 (a special case of ®) such that its action between / and 
an axisymmetric function h (satisfying 7 ^(o,o, 7 )^ = V 7 ) yields a function on S^x M+ defined by 

{fQh){n',t) = ifin^^Tth) = [ dn{n) f drr^ f{n,r) {TZ^^Tt gnn,r), (46) 

Js2 Jr+ 

where Tin' is the rotation operator for axisymmetric functions, i.e., Tin' = ,9',o)- what follows. Tip and Tin refer 

to directional and axisymmetric rotations, respectively, since p and n denote angles on SO(3) and §^, respectively. 
More detailed discussions about rotation, translation, and convolution operators on the sphere and the radial line 
can be found in Refs. [46] and [39, 68 ], respectively. All operators can be evaluated exactly via their harmonic 
representations and all transforms from pixel to harmonic space and conversely are exact thanks to the sampling 
theorems adopted. 


IV. SPIN FLAGLET TRANSFORM 

We now extend the scalar, axisymmetric flaglet transform of Ref. [39] to analyze signals of arbitrary spin and to 
probe their features in the angular direction. This new separable flaglet transform is defined in Sec. IV A, where we 
present the core equations and rely on the previous translation and convolution operators. However, to highlight the 
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generic properties of the transform we do not specify the details of the flaglet kernels and until Sec. IVB, 
where we also present our efficient implementation using sampling theorems and multiresolution algorithms. 


A. Transform 


The flaglet coefficients of a spin signal sf ^ x M+), denoted by \ are defined on SO(3)x]R+ as the 

directional convolution of g/ with the flaglets G x M+) 

wf'{p,t) = U®s^'^){p,t) = [ dll(n) / drr\Jin,r) (UpTi (47) 

where i = Iq, ... and j = Jq, ... denote the angular and radial scales, respectively. Flaglets are localized in 
scale, position and angular orientation, and designed to capture the directional, high-frequency content of the signal 
of interest. A scaling function G x IR+) is introduced to capture the low frequency content of the signal in 

the scaling coefficients . In this work we analyze the low frequency part of the signal using axisymmetric flaglets 
because its directional content is typically of lower interest. In the context of weak-lensing, largest scales (low k in 
Cf^) have far fewer modes, contain less information (due to the reduced amplitude of Cf^)^ and are more difficult to 
deal with due to observational masks and other considerations mentioned above. However, extending the formalism 
presented here to support directional scaling functions is possible and straightforward. 

We define G I/^(§^x M+) as the axisymmetric convolution of g/ with the axisymmetric scaling function 

W^f{n',t) = = [ dfi(n) [ drr^ ,/(n,r) ,$)*(n,r). (48) 

Since is axisymmetric, the Fourier-Laguerre coefficients s^£m,p are nonzero for m = 0 only. 

In what follows we will use a compressed notation for all sums over scales i and j, assumed to run over Iq, ... 
and Jo,..., J, respectively. These bounds, as well as the flaglets and scaling functions and are consistently 
defined in the next section in order to achieve exact reconstruction of the input (bandlimited) signal. 

The wavelet transform makes direct use of the rotation, translation and convolution operators defined previously. 
Hence, the forward flaglet transform (or analysis step, equivalent to Eq. (47)) can also be expressed in Wigner-Laguerre 
space (on SO(3)xM+) in a straightforward manner by 


(w:f) = 

\ / mn,p 

The scaling coefficients (on x IR+) in Fourier-Laguerre space read 


8^^ . ^ij * 

21 + -^sj£m,p s^£n,p 


s{W:f)em,p = {W:fUYe^Kp) = ^ 


s^£m^p • 


(49) 


(50) 


The reconstruction of the original signal is achieved through the inverse flaglet transform (synthesis step) 

sf{n,r)= f dn{n')[ drr^ W‘f{n') s<^){n,t) + V / dp{p) [ drr^ W‘f" (p) {11^% (51) 

Jso{3) Jr+ 

or in Fourier-Laguerre space 

siW:f)e^,p + EE (^2) 


ij n 


In order for the transform to be exact, i.e., for the wavelet coefficients to capture all the information content of sf, 
the flaglets and scaling function must be defined to satisfy Eqs. (49), (50), and (52). In other words, they must be 
chosen such that their Eourier-Laguerre coefficients satisfy the admissibility condition 


Att 

2 ZT 1 


\s^m,p\^ 


87r2 

2i+l 


Ei^^Vpi' = i- 


\/£,p . 


(53) 
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Tiling of the angular harmonic space i Tiling of the radial harmonic space p 



FIG. 2 : Tiling of the angular (left panel) and radial (right panel) harmonic lines employed to construct flaglets in Eq. (59). 
The left panel shows the kernels /^a(^) and tiling the range ^ = 0,..., L — 1 (zero is not shown but is included in px), 

while the right panel shows Ki, and pi, tiling p = 0,..., P — 1. Since L = P = 243 and A = z/ = 3, six kernels are needed to fully 
tile [0, L — 1] and [0, P — 1], i.e., the scales considered are i = 0,..., 5 and j = 0,. .., 5. The first two scales are incorporated in 
the scaling function by setting /q = Jo = 2. The thicker lines highlight the scales shown in Figs. 3 and 4. Note that we do not 
show the hybrid kernel pxu, included in the scaling function to satisfy the admissibility condition at low p and i, see Eqs. (53) 
and (61). We also do not include the directional component Qm in this figure, which is nevertheless shown in real space in 
Fig. 3. 


B. Wavelet construction 


So far we have constructed a generic flaglet transform relying on the properties of the Fonrier-Lagnerre space. To 
uniquely characterise this transform, we need to specify the scaling function the flaglets and the bounds 

Jo, /, Jq: J to satisfy Eq. (53). We follow the construction of scalar flaglets [39] and scale-discretized wavelets 
[46, 70]. Other types of wavelets could be used; we opted for scale-discretized wavelets because they exhibit good 
localization properties [71] and can be easily extended to probe directional features [46, 70]. They exhibit a compact 
representation in harmonic space that also allows several computational improvements, which are essential to make a 
3D spin directional transform tractable. We use tiling functions n and p defined as 


where 


Kx{t) = x/k\{t/\) - kx{t), 

= Vkx{t), 


kx{t) 




(54) 

(55) 


(56) 


which is unity for t < 1/A, zero for t > 1, and is smoothly decreasing from unity to zero for t G [1/A, 1]. In these 
equations, A is the wavelet dilation parameter and characterises the density of the tiling. In what follows, A will refer 
to the wavelet tiling in the angular direction, while a second parameter u is introduced to tile the radial direction 
(with functions and p^, constructed in the same manner). Finally, the remaining free function sx{t) (also defined 
for u) must have compact support in [y, 1] and is defined as — 1/A) — 1) with 


s{t) 


e t e [-1,1] 

0 , tii-1,1] ■ 


(57) 


The tiling of the radial and angular harmonic spaces constructed with these functions is shown in Fig. 2 for A = = 

3, Jo = Jo = 2, L = P = 243. 

In addition, as for the scalar axisymmetric flaglet transform [39], a hybrid tiling function is needed in order to 
construct a suitable scaling function <F and satisfy Eq. (53), 


t')=x/kx{t/\)K{t') + kx{t)K{t'/u) - kx{t)K{t'). 


( 58 ) 
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With the tiling functions n and r] defined in Eq. (54) and Eq. (55), we construct the flaglets in Eourier-Laguerre 
space as 


Om. (59) 

The extra function ( controls the directionality component [47, 48, 70] and is parameterised by an azimuthal band 
limit N, such that Qrn = 0, V^, m with \m\ > N. It is defined in harmonic space as 


Cim — 




1 

¥ 




(60) 


to satisfy various azimuthal symmetries and normalized such that ¥m ~ values of £ for which (irn are 

nonzero for at least one value of m (see, e.g., , [71]). In this expression, a = 1 for — 1 odd and a = i for — 1 even, 
6 = [1 - {-l)^+^]/2 and c = min{A^ - 1,^ - [1 - (-1)^+V2}. 

Einally, the axisymmetric scaling function is needed to capture the information from the regions £ < and 
p < ^ which are not probed by the flaglets defined above. Thus, it is constructed to satisfy Eq. (53) and reads 


= 




0 . 


U-^O ) ’ 


if ^ > A^o, p < 

\i£< A^o, p > 

\i £ < A^°, p < 

elsewhere. 


(61) 


These definitions entirely characterise the spin, directional flaglet transform. The free parameters of the transform 
are A, z^, A^, Jq, and Jq. A and v define the support of the flaglets in Eourier-Laguerre space and therefore the scales 
and features the flaglet coefficients extract, as shown in Eig. 2. I and J are the maximum angular and radial fiaglet 
scales captured, respectively, and are fixed by the radial and angular band limits of the signal: / = [log;)^(I/ — 1)] and 
J = [log^(P — 1)], to satisfy Eq. (53). N is used to specify the number of directions to probe. Jq and Jq are the first 
angular and radial scales of interest (larger scales being captured by the scaling function) and can be freely chosen 
provided they satisfy 0 < Jq < / and 0 < Jo < J. 

The flaglets (and scaling function) tile the angular and radial frequency domains £ and p, as illustrated in Eig. 2, 
while satisfying the admissibility condition of Eq. (53). Wavelets are well localized simultaneously in the spatial 
domain, both in position and orientation, and the harmonic domain. It is shown in Ref. [71] that scale-discretized 
wavelets on the sphere exhibit excellent concentration properties, both in the scalar setting and in the spin setting. 
They are also steerable [47, 48, 70]. Elaglets naturally inherit these properties. In particular, as shown in Ref. [68] for 
the scalar setting, the flaglet transform forms a tight Parseval frame, i.e., the norm of the input signal is conserved. 

Eigs. 3 and 4 show spin-0 and spin-2 flaglets resulting from our construction, with parameters X = u = 3 and 
^0 = 'A) = 2. In contrast to the Eourier-Laguerre basis functions shown in Eig. 1, which are not localized in real 
space (but are delta functions in harmonic space), all flaglets have good spatial localization properties in both radial 
and angular dimensions. They are also localized in Eourier-Laguerre space by construction, as shown in Eig. 2 and 
Eq. (59). Note that the particular i = 2, 3 and j = 2, 3 scales shown in Eigs. 3 and 4 are highlighted as thicker lines in 
Eig. 2. Eig. 3 shows spin-0 flaglets (real part only), comparing the N = 2^3 flaglets to the axisymmetric N = 1 case 
(first shown in Ref. [39]). The angular components have even and odd symmetry for N = 2,3, respectively, by the 
construction of ( [46]. Eig. 4 shows the real and imaginary parts and the complex modulus of spin-2 axisymmetric 
{N = 1) flaglets. 


C. Multiresolution algorithm and implementation 

As described previously, all the integrals and convolutions related to the Eourier-Laguerre and Wigner-Laguerre 
transforms can be computed exactly via their harmonic space representations and by appealing to sampling theorems 
to analyze or reconstruct signals in pixel space. Let us consider an input signal sf ^ I/^(S^x IR+) with radial and 
angular band limits P and L, respectively. As before, such a signal can be represented exactly by ^ 2PL‘^ samples 
and PL^ spherical harmonic coefficients. The full algorithm to perform the flaglet transform is described by the 
pipeline below. The first and the third operations are spin Eourier-Laguerre transforms, while the second is the 
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FIG. 3: Real part of spin-0 directional flaglets, plotted up to r = R/2 (with R the last node of the radial sampling scheme), 
showing their excellent localization and directional properties (even and odd for iV = 2, 3 respectively). This demonstrates 
that the flaglet transform can probe radial and angular modes well defined in both real and harmonic space. The radial and 
angular harmonic modes i and p corresponding to the scales i = j = 2,3 are highlighted as thicker lines in Fig. 2. The dashed 
lines and half sphere show the slices that are plotted in the various subpanels of this figure. Flaglets are localized in both pixel 
and frequency space, unlike the Fourier-Laguerre basis functions shown in Fig. 1 which are delta functions in harmonic space 
and therefore not localized in pixel space. 


wavelet transform harmonic space filtering operations of Eqs. (49), (50), and (52). 


X Fourier-Laguerre « Wavelet 


I (W‘f) , (wf)' I 

\ sj J \ J rnn,p J 


Fourier-Laguerre 
^^^ 




The flaglet transform will produce flaglet coefficients on SO(3) xM+ for (/ — Jq + 1)( J — Jo + 1) scales, each captured 
in ^ 4NPL^ samples for an azimuthal band limit N (controlling the number of directions). The additional scaling 
function is axisymmetric and captured in ^ 2PL^ samples on §^x M+. Therefore, the total number of coefficient scales 
is Jtot = (4 — /o + l)(J — Jo + l) + l- Each one requires the filtering of Eq. (49), scaling as 0{LP‘PN), and a Wigner- 
Laguerre transform (or, for the scaling coefficients, filtering of Eq. (50) and a Eourier-Laguerre transform). Therefore 
the Wigner-Laguerre transforms on SO(3) xM+ will dominate the computation time, and the overall complexity of the 
flaglet transform is 0{JtotNPL^). Note that for N > 1 the scaling function computations are negligible compared to 
the wavelet computations, in which case we have Jtot = (4 — /q + 1)(J — Jo + !)• 

However, by definition of the flaglets and transform (Eqs. (59) and (52)), the flaglet coefficients have lower band 
limits than the original signal. Specifically, 


^ 0 and (Wpy ^ 0 only if £ G [V-\ m = n =-N,... ,N, p G [,,■ 

' 777 / 77 / 


J-l ;J + 1 


p2) 


In other words, the flaglet coefficients ^ have angular and radial band limits and respectively. Thus, 

the number of samples needed to capture the information is much smaller than at the full resolution. One can employ 
a multiresolution algorithm and use the minimum band limits and number of samples for each flaglet coefficient 
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FIG. 4: Same as Fig. 3, but for spin-2 flaglets and showing their real and imaginary parts, as well as the complex modulus. 
This demonstrates that the flaglets natively support the symmetries of spin signals while being well localized in both real and 
harmonic space. 


scale (and similarly for the scaling coefficients). This reduces the complexity of the fiaglet transform significantly, 
to OiNPL^), since the computation of the largest scales ij now dominates that of all other scales and the scaling 
function. 

We have implemented the exact and efficient algorithm described previously to compute the spin directional fiaglet 
transform in the existing FLAG and FLAGLET codes (http://www.flaglets.org). The latter relies on the SSHT 
(http://www.spinsht.org), s2let (http://www.s21et.org) and S03 (http://www.sothree.org) codes for the 
angular transforms and sampling theorems. The EFTW (http://www.fftw.org) code is used to compute Fourier 
transforms. The core algorithms of the Fourier-Laguerre, Wigner-Laguerre and fiaglet transforms are implemented in 
C and are able to handle large band limits and billions of samples on IR+ and S0(3)x]R+. Interfaces to the core 
functions as well as convenient data manipulation and plotting routines are provided in Matlab and Python. 


V. APPLICATION TO 3D WEAK TENSING 

In previous sections we developed the spin fiaglet formalism to extract content localized in both space and frequency 
from 3D signals of arbitrary spin, which has the added benefit that radial and angular modes are separable. This 
is a general framework that can be applied to arbitrary 3D spin signals. We now apply this approach to represent 
3D cosmic shear and other weak leasing observables, and relate the covariance of their fiaglet decompositions to the 
leasing power spectrum. In particular, we highlight how one can exploit the properties of flaglets (separability and 
simultaneous localization in frequency and pixel space) to deal with partial sky coverage and small-scale modeling 
uncertainties. We also examine the approximations commonly adopted in cosmic shear analyzes and apply them to 
the fiaglet covariance of weak leasing observables. Note that we focus on using axisymmetric wavelets to increase 
simplicity and readability of the final estimator. However, the latter is straightforwardly extended to directional 
flaglets, and this does not affect any of the advantages and properties discussed below. 
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A. 3D cosmic shear 

We relate the covariance of the flaglet representation of the cosmic shear field to the leasing power spectrum. As 
shown in Sec. II, the 3D cosmic shear signal 27 is characterised by the leasing power spectrum k') defined in 

Fourier-Bessel space by 


( 2 lem{k) (63) 

where all relevant physical effects are assumed to be modelled (as discussed in Sec. II C). As in Sec. II, the brackets in 
this section denote ensemble average over cosmological realisations, rather than inner products (which are distinguished 
by a vertical bar between the functions of the inner product). 

We now consider the wavelet transform of 27 , using axisymmetric spin flaglets {i.e., N = 1). In this case all 
directional convolutions ® become axisymmetric convolutions ©. The flaglet coefficients live on x IR+ and read 


W^^'\n,r) = ( 27 © 2 ^'*-^)(n,r). 
We consider the covariance between flaglet scales. 


^ (n, n', r, r') = ( (n, r) 




(64) 


(65) 


Injecting the Fourier-Bessel decomposition of 27 into the integral form of the convolution of (n, r), one obtains 

^ [ dfcfcT dk'k'^Cf‘^ik,k') 2F^iik,n,r) (66) 


where wavelet transform of the Fourier-Bessel basis functions: 

2F7(fc,n,r) = {2Ye^Mk-)Q2^'^){n,r) (67) 

= [ d^in') [ drY'\YUn')jeikr') {nr,Tr2^'^)^n',r'). (68) 

J^+ 

Note that this formalism can be straightforwardly extended to directional wavelets: directional convolutions then 
replace axisymmetric convolutions, and the angles n, n' G become p, p' G SO(3) since the flaglet coefficients then 
live on SO(3)x]R+. 

The previous expression can be simplified by appealing to the harmonic representation of the translation and 
rotation operators, which when applied to axisymmetric flaglets 2 ^*-^ (^-^v Refs. [46, 72]) yield 


47r 


(TrlZr^ ^ ^^K;ir)Y;M 


The functions 2 F^^{k,n,r) then take the simple form 


2 FZ{k,n,r) = 2 %^{k,r), 


(69) 


(70) 


where 


2^77 ik, r) = J2 JeAk) KAr) 2^70,p- ( 71 ) 

P 

As shown previously, the functions admit an analytical form [39] and can be tabulated. Furthermore, 2 ^Wk^r)^ 
and consequently 2 Tg*^(/c, n, r), is straightforward to compute since only a one dimensional summation over a finite 
number of p samples is required. 
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Further simplifications can be made by exploiting the spherical harmonic addition theorem: 

YeMYl^in') = ■ n'), (72) 

m 

where Pi are the Legendre polynomials. By the additional theorem one can infer that the wavelet covariance only 
depends on the angle A6>, where n • n' = cos(A6>). Thus, ^ (n, n', r, r') becomes 

f dk'k'^ Cf^ik,k') Pe{A0) 2 nyk,r) (73) 

^ ^ 4 J^+ J^+ 

The angular dependence of the flaglet covariance depends only on the angular separation between pixels on the sky. 
This is expected since we are analyzing the 2-point fluctuations of the signal with axisymmetric flaglets and have 
assumed statistical homogeneity and isotropy. 

Eq. (73) is one of the key results of this paper. One can build a full likelihood analysis of 3D weak leasing observables 
and constrain cosmological models and parameters using spin flaglets (similar to that of Ref. [12]). On the left hand 
side of Eq. (73) is a quantity that can be computed from data using a wavelet decomposition (here using axisymmetric 
flaglets, but a similar expression can be obtained with directional flaglets). On the right hand side is the theoretical 
expectation of this quantity related to cosmological parameters through the 3D lensing power spectrum 
Even though Cf^{k^k') mixes radial and angular modes, these are isolated into scales by the kernels 

2 'HWk^r). In practice, the flaglet covariance can be calculated for a pair of flaglet coefficients ij and i'j' and a pair 
of radii r and r' by averaging over pairs of pixels separated by AO. This approach naturally accounts for the 2+ID 
nature of the shear observables and takes advantage of the convenient connection between the Eourier-Laguerre and 
Eourier-Bessel transforms. It takes full advantage of the separability of the Eourier-Laguerre and flaglet transforms. 
Going through flaglet space allows one to exploit the excellent localization properties of flaglets in both pixel and 
frequency space, a property absent in all alternative approaches. Simultaneous localization in both pixel and frequency 
space allows one to cut or filter regions of the sky due to unobserved, unreliable or contaminated data, at the same 
time as cutting or filtering harmonic modes. Eor example, one may wish to cut or filter small-scale harmonic modes 
where physical modeling is less certain, e.g., high-/c (and k') modes from Cf^{k^ k'). This is possible by studying the 
kernels 2 kCi {k^ r) and only considering the flaglet scales ij, i'j' that probe the physical scales of interest in Cf^{k^ k'). 
Thus, the spin flaglet 3D weak lensing approach is a flexible way to avoid the complications of correlation functions 
and Eourier-Bessel representations, which are typical of standard approaches. 


B. Generalization to other weak lensing quantities 


As shown in Eq. (6), spin-2 shear distortions are not the only effect of matter fluctuations on the observed properties 
of galaxies [49]. Other observables are expected to be measured at much lower signal-to-noise ratio than the shear 
distortion but they nevertheless contain important cosmological information. Thus, we generalize the flaglet covariance 
derived previously to other lensing distortions of different spin. Eor this purpose we introduce a general scaling factor 


-ivW+iF, 

\Wt{£ + 1){£ - me + 2)2 + 2yp{£ + 1 ) 3 ] 


Mi^s 


1 — I AT 

2VC^ “ 


(^+3)! 


= -hNi;^ 


if s = 0 
if 5 = 1 

if 5 = 2 
if 5 = 3 


2Y (^-3)! 2^ 

SO that we can write sXim{k) = Mi^s4^im{k), with defined in Eq. (6). In this unified approach, we have 


(n, n', r, r') = — 

TT 


(74) 


Y,iMe,sy [ dkk^f dk'k’yf^{k,k') yy{k,n,r) y2*(k',n'y) (75) 
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with sF^^{k,n,r) generalizing Eq. (70) to other values of spin. As before, the extension to directional wavelets is 
straightforward. Exploiting the spherical harmonic addition property, again, leads to the simplification 

Cw'i'(A6i,r,r') = 1 V [ dkk^ [ dk'k'"^ Cf‘^{k,k') Pe{Ae) {k,r) sn^^'*{k',r'). (76) 

^ ^ ’ Jr+ Jr+ 

with s'kiWk^r) defined as previously. This expression generalizes the 3D cosmic shear power spectrum to the case 
of lensing distortions of other spin. As before, this quantity can be calculated from data by averaging over pairs of 
pixels separated by AO in bins of ijf, i'j', r, r'. The theoretical prediction through the lensing potential, the Legendre 
polynomials and the s'kL kernels is readily computable. Einally, this formalism can also be extended to support 
cross-terms between the various spin components, which would be required for a combined analysis. 


C. Approximations 

We have related the covariance of weak lensing observables in fiaglet space to the underlying 3D lensing power spec¬ 
trum in the general setting. Several approximations are used in the literature to make the computation of shear and 
other lensing quantities more approachable. Some of these are lossless—or nearly so—while others result in a loss of 
information with respect to the general case. The most commonly used approximation is to use cosmic shear ‘tomog¬ 
raphy’, in which flat-sky and Limber approximations are typically also made. We apply these approximations below 
to the fiaglet covariance of spin lensing quantities, where appropriate, and show that they result in representations 
that are more readily computable but are a lossy representation of the 3D lensing power spectrum. 


1. Flat-sky approximation 

Weak gravitational lensing of galaxy images has the strongest signal-to-noise on the scales of groups or clusters 
{e.g.y Ref. [55]) and therefore there is relatively little signal in the power spectrum on large angular scales. In 
addition, for surveys with relatively small sky coverage and/or surveys with several small observational fields, effects 
due to the spherical geometry of the setting are expected to have a relatively small impact. However, for future 
surveys {e.g., Euclid [73, 74] and LSST [75]) with observations over increasingly greater coverages of the sky, flat-sky 
approximations will become increasingly less accurate. 

The flat-sky approximation assumes that the angular extent of the observational field is small and hence the 
geometry of the angular component is assumed to be planar (z.e., Euclidean). In this setting the spherical harmonics 
are approximated by the product of Bessel functions (of the first kind) and complex exponentials [7, 76]. In the 
case of Eqs. (73) and (76), however, the sum over the product of spherical harmonic functions reduces to a Legendre 
polynomial (due to the spherical harmonic addition theorem), which in the fiat-sky case would be approximated 
by a zeroth order Bessel function. However, given the computational ease of computing Legendre polynomials the 
full-sky case is readily computable and so a flat-sky approximation is not required. This differs to the computation of 
Cf^{k^ k')^ where a fiat-sky approximation can reduce computational time considerably and has so-far been assumed 
in current applications to existing data [11, 12], where sky coverage is relatively small. 


2. Limber approximation 

The Limber approximation [77] assumes that the evolution of radial modes is small over the survey volume under 
consideration. Under this approximation the spherical Bessel functions may be approximated by Dirac delta functions: 

k^/‘^ji{kr) - + 1/2 - kr), (77) 

as shown in Refs. [8, 59], and rederived in Appendix B for the spherical Bessel convention adopted herein. 

The use of this approximation simplifies matters considerably: an 8-nested integral equation for Cf^{k^ k') becomes 
a single integral (see, e.g., Refs. [9, 12]). However, even though this approximation primarily affects large scales 
£ < 100, it is severe and can cause errors in the inferred cosmological parameters of tens of percent [9]. Exploiting 
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the Limber approximation simplifies the Ji^p{k) functions to 


Je,p{k) 


(£+1/2)3/2 .^+i/2x 

\j 2 fc7/2 n k )' 


( 78 ) 


and the 2 'Hl^{k,r) kernels to 


2 'Hy {k, r) 


(£+1/2)3/2 
V 2 fc7/2 




(79) 


This is a significant simplification since, although Ji^p{k) already does not need to be computed by the integration of 
products of Bessel and Laguerre functions (see Eq. (40)) and can be computed analytically (as shown in Ref. [39]), 
the analytical computation can also be challenging for high-^. Eq. (78) and, consequently, Eq. (79) can be readily 
computed both accurately and efficiently for all 


3. Tomographic approximation 

The tomographic approximation is a combination of the the fiat-sky and Limber approximation with the addition 
of a discretisation of the signal along the radial direction into a series of redshift bins. Within each redshift bin the 
projected 3D contribution, from a range of galaxies with redshifts assigned to that bin, is used to construct a 2D 
power spectrum. This binning in redshift is lossy as any evolution that occurs in the signal on scales smaller than the 
bin widths is lost. Nevertheless, due to computational and conceptual ease it is the most widely used approach. 

To relate the wavelet decomposition of the data to the tomographic power spectrum we do not require radial 
transforms but can simply relate the wavelet covariance to the cosmic shear tomography power spectrum using the 
approximation 


Cf‘l’{k, k') ~ + 1/2 - kr)5°{£ + 1/2 - fcV'). 

By this tomographic approximation, Eq. (76) reduces to 

^4 


(80) 


C-7-'/(i»,.,/) criry) P,m (81) 


where the kernels are those given in the Limber-approximated case of Eq. (79). The final approximation for 
tomography with binning in redshift is 


{A0,rhi„{za),rhin{za')) ^ [ dr [ dr'{A0,r,r') W{r,r',r],in{za),r],in{za')), (82) 

Jr+ Jr+ 

where the weight function W describes if a galaxy is in the bin-pair combination (rbin(^a )5 ^bin(^a')) 
weight functions are usually designed as top-hat functions in spectroscopic redshift Za and corrections are applied for 
‘leakage’ due to photometric redshifts (see, e.g., Ref. [9]). The methods presented here can also be extended to deal 
with redshift uncertainties and more optimal weight functions. 


VI. CONCLUSIONS 

We have developed a novel approach to analyze weak gravitational lensing observations using 3D wavelets (flaglets) 
and the Eourier-Laguerre transform. We have shown that the covariance of the flaglet coefficients of weak lensing 
observables (such as cosmic shear, size distortions and flexion) can be directly related to the 3D lensing power spectra, 
exploiting the analytical connection between the Eourier-Bessel and Eourier-Laguerre transforms. Thanks to the 
separability of their radial and angular components, the Eourier-Laguerre and flaglet transforms can take advantage 
of existing fast algorithms and sampling theorems developed for the radial line, the sphere, and the rotation group. 
Eurthermore, this approach is suited to dealing with the 2+ID nature of weak lensing observations and cosmological 
data sets in general. Eor example, typical data complications and systematic uncertainties are given on the sky and in 
the radial direction separately. The excellent localization properties of flaglets in both pixel and frequency space can 
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be exploited to deal with these effects and alleviate the complications of weak leasing analyzes {e.g., complicated sky 
masks and uncertainties of the small-scale modeling) without complicating the estimation of the ffaglet coefficients 
and covariance. Therefore, the method introduced here is a powerful alternative to common real and harmonic-space 
approaches where these effects are difficult to deal with [11, 12]. Leasing observables and covariances are more easily 
modelled in the Fourier-Bessel basis, where unreliable or unmodelled small scales can also be filtered out. However, 
the mode-mixing due to partial sky coverage is difficult to handle in Fourier-Bessel space. By contrast, real space 
correlation function measurements can naturally deal with complicated survey window functions, but filtering scales 
and modeling covariances is then significantly more difficult. The dual localization of wavelets in pixel and harmonic 
space gives them the advantages of both approaches. In future work we will demonstrate the ability of this new 
approach to robustly extract cosmological information by applying it to simulations and real observations of weak 
leasing signals. 

We have updated the publicly available FLAG and flaglet (http://www.flaglets.org) codes to compute the 
spin Fourier-Laguerre and flaglet transforms. These rely on the fast algorithms implemented in the following codes, 
also publicly available: s2let (http://www.s21et.org) for the spin directional spherical wavelets, SSHT (http:// 
WWW. spinsht. org) for the spherical harmonics transform, S03 (http: //www. sothree. org) for the Wigner transform, 
and EETW (http://www.fftw.org) for the Fourier transform. In all these codes, the core routines are implemented in 
C, and we have provided numerous wrappers and interfaces in Python, Idl, and Matlab, to call the core routines, 
and to manipulate and visualize data sets. As highlighted in Ref. [39] for the scalar setting, these codes can deal with 
large band limits and large 2D and 3D data sets, which can be pixelized into billions of samples and manipulated 
with no loss of information. 
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Appendix A: Properties of special functions 

We concisely review the orthogonality and completeness properties of the special functions considered throughout 
this paper, of which we make continual use. In particular, we consider the spherical Bessel functions, the Laguerre 
polynomials, the spin spherical harmonics, for which we also review the spin raising and lowering operators, and 
finally the Wigner functions. 

The spherical Bessel function (of the first kind) of order denoted by is defined by 

k{x) = (Al) 

where is the standard Bessel function of the first kind. The closure relation for the spherical Bessel functions is 
given by 

J^ drr‘^je{kr)je{k'r) = - k'), (A2) 

which plays the role of both completeness and orthogonality relations in the continuous setting (note that the spherical 
Bessel transform does not admit a sampling theorem leading to exact quadrature [39, 68]). The spherical Laguerre 
basis functions Kp are related to the standard Laguerre polynomials by Eq. (24). The orthogonality of the spherical 
Laguerre functions reads 

{Kp\Kg)= [ drr^Kpir)K;ir) = 5f^, (A3) 

Jr+ 

while completeness is inherited from the completeness of polynomials [39]. Both the spherical Bessel and Laguerre 
functions form a complete basis for the representation of functions defined on the radial line M+. 
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The orthogonality and completeness of the spin spherical harmonic functions sYim read 

isYemUYem') = [ dn(n),YUn)sY;,^,(n) = 


(A4) 


and 


'^sYem{n)sYl^{n') =S^{n-n') = S^{cose - cose')S^{(p - (j)'), (A5) 


respectively. The spin spherical harmonics are the canonical (complete) basis for the representation of functions 
defined on the sphere while accounting for the symmetry of spin-5 signals. Spin raising and lowering operators, 
5 and 8 respectively, increment and decrement the spin order of a spin-5 function and read [51-54] 


= — sin^ 0 


d 


i d 


dO sin 0 i 


-sin ^ 


d 


i d 


dO sin 0 d(j) 


sin-" 6>, 
sin" 0, 


(A6) 

(A7) 


respectively. Using these operators, spin-s spherical harmonics sYtm can be constructed from scalar (spin-0) harmonics 
by 


sYim{n) = _s d^Yimin), for 0 < s < £ (A8) 

sYim{n) = Ni^+s d~‘'Yim{n), for -£<s<0. 

We assumed the Condon-Shortley phase convention and define the factor 


Ne,s 


(^ + s)! 

i£-s)\- 


(A9) 


The orthogonality and completeness of the Wigner functions read 

= [ d^i{p)DlMDiln'iP) = (AlO) 

JSO(3) + i 

and 

“ P') = “ 0 ;')^^ (cOS/3 - COS/3')S^(j - j'), (All) 

£mn 

respectively. The Wigner functions provide a complete basis for the representation of functions defined on the rotation 
group SO(3). Note that we adopt basis functions since this simplifies connections to wavelet transforms. 


Appendix B: Limber approximation 


The extended Limber approximation is derived in Ref. [77]. By comparing Eqs. (5) and (12) of Ref. [77] and 
neglecting higher order terms, we find: 

J^ dkk‘^ji{kr)ji{kr')F{k) ~ (Bl) 

which was noted previously by [8]. It follows that the Limber approximation can be represented by the following 
approximation of the spherical Bessel functions: 




5^{£+ 1/2-kr). 


(B2) 
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We note a difference of a factor of compared to [8] and an additive factor of 1/2 in the argument of F in Eq. (Bl). 
Note that [8] consider a different convention for the spherical Bessel transform to that considered herein (which differ 
by a factor of F). We follow the same convention as [6]. 
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